Dynamics of confined water inside carbon nanotubes based on studying tetrahedral order parameters

Water dynamics inside hydrophobic confinement, such as carbon nanotubes (CNTs), has garnered significant attention, focusing on water diffusion. However, a crucial aspect remains unexplored - the influence of confinement size on water ordering and intrinsic hydrogen bond dynamics. To address this gap, we conducted extensive molecular dynamics simulations to investigate local ordering and intrinsic hydrogen bond dynamics of water molecules within CNTs of various sizes (length:20 nm, diameters: 1.0 nm to 5.0 nm) over a wide range of temperatures (260K, 280K, 300K, and 320K). A striking observation emerged: in smaller CNTs, water molecules adopt an icy structure near tube walls while maintaining liquid state towards the center. Notably, water behavior within a 2.0 nm CNT stands out as an anomaly, distinct from other CNT sizes considered in this study. This anomaly was explained through the formation of water layers inside CNTs. The hydrogen bond correlation function of water within CNTs decayed more slowly than bulk water, with an increasing rate as CNT diameter increased. In smaller CNTs, water molecules hold onto their hydrogen bond longer than larger ones. Interestingly, in larger CNTs, the innermost layer’s hydrogen bond lasts a shorter time compared to the other layers, and this changes with temperature.


Methods
MD simulations were performed under a constant number of particles, pressure, and temperature ensemble to explore the ordering and the HB dynamics of water confined within CNTs of varying sizes.Water molecules were confined within arm-chair CNTs with diameters of 1.0 nm, 1.4 nm, 2.0 nm, 3.0 nm, and 5.0 nm, respectively.The chirality indices correspond to diameters are (8,8), (10,10), (15,15), (22,22), and (37,37), respectively.The initial configuration of the water-CNT system is depicted in Fig. 1.Notably, the length of the CNT was set to 20 nm for these simulations.

Simulation setup
MD simulations were performed using NAMD package 34 , employing the SPC/E water model 35 , known for accurately predicting various bulk water properties.Non-bonded interactions between carbon atoms were described using the Lennard-Jones potential with parameters from Werder et al. 36 .CNT remained fixed throughout the simulations, imposing a positional restraint on each carbon atom.The system was immersed in a periodic simulation box of 4 × 4 × 23 nm 3 .The dimensions of the simulation box vary along X and Y directions with the CNT diameter.The confined water inside the CNT was investigated in the temperature range of 260 K to 320 K, with increments of 20 K.The Langevin thermostat set the temperature to the desired value and a Nose-Hoover Langevin piston maintained the pressure at 1 bar, with a period of 100 fs and a damping time scale of 50 fs.Unlike previous studies with infinitely long CNTs along the axis, our simulation allowed water molecules Figure 1.Schematic representation of Water-CNT used in the MD simulations: (A) 1.0 nm (8, 8), (B) 2.0 nm (15,15), (C) 3.0 nm (22,22), and (D) 5.0 nm (37, 37) of CNT filled with the water molecules.Oxygen, Hydrogen, and Carbon atoms are shown in red, white, and green colors.The chirality index of each arm-chair CNT diameter is given in brackets.
to flow in and out of CNTs with a finite length of 20 nm, which gives a closer representative of real samples of CNTs used in various experiments.The Particle Mesh Ewald summation method (PME) was used to compute long-range electrostatic interactions.A time step of 2 fs is used for the integration of the equation of motion.Bonded interactions were computed at every time step, while non-bonded interaction was calculated every two steps, with a cutoff of 12 Å and a switching function of 10 Å.All the simulated systems were minimized for 10000 steps, followed by equilibration at the target temperature for 50,000 steps (100 ps) before the 50 ns production run.The system configuration was saved every 500 steps (1 ps) for subsequent analysis.

Local order parameters for water
To characterize the arrangement of water molecules confined inside various sizes of CNTs, we have computed the order parameters that describe the water structure locally.Typically, confined water structures are mostly investigated using molecules density profile.In this work, we also used order parameters: Orientational Tetrahedral Order ( S q ) and Transnational Tetrahedral Order ( S k ).These parameters offer a distinct advantage over density profiles as they provide deeper insights into the structural order of water within the hydrophobic confinements of CNTs.The S q order parameter is given as: ψ ij , is the angle between lines joining the oxygen atom of the water molecule under consideration and its nearest neighbor oxygen atoms i and j.The angle between one atom to the four nearest neighboring atoms is measured and averaged.S q equals to 1 when cosψ ij =-1 3 , which is a regular tetrahedron.S q is widely used to describe the structure of liquid water, measuring the degree of flexibility in the arrangement of water molecules.S q value of 0 corresponds to an ideal gas, while a value of 1 corresponds to the regular tetrahedral structure of ice.
The S k order parameter is described as: r k is the distance between the oxygen atom in water molecules and its nearest neighbor of oxygen atoms.r is the average distance with the four nearest neighbor oxygen atoms.S k relies on the variance of radial distances between the oxygen atom of a water molecule under consideration with four nearest oxygen atoms.S k equals 1 corresponds to the regular tetrahedron.Water molecules near the CNT walls are excluded from the order parameter calculation due to their proximity to carbon atoms, which prevents them from occupying the center of a tetrahedral configuration.order parameters were calculated over a whole trajectory length i.e.50 ns.Error bars represented by the standard deviations were estimated by calculating the order parameter for every 5.0 ns at a sampling rate of 1.0 ps.

Hydrogen bond correlation function
The dynamic and thermodynamic characteristics of water are intricately connected to the nature of hydrogen bonding.Therefore, it is essential to gain a thorough understanding of the behavior of HBs in water molecules at various temperatures and within CNTs of different sizes.In order to investigate the dynamic processes involved in forming and breaking the HB network in a confined environment, we have calculated the HB auto-correlation function (HBACF) using the following geometrical criterion: where b(τ ) is the HB order parameter.For b(τ ) , a value of 1 is assigned if the two assigned water molecules form an HB at time τ , and 0 if no HB is formed.The decay of HBACF is relatively swift, typically occurring within a few picoseconds, much faster compared to other auto-correlation functions for other simple liquids.The HBACF provides a measure of the likelihood that a pair of water molecules, initially bonded at time τ , will remain bonded at time t + τ .In other words, the HBACF provides insight into how persistent and dynamic the HB network in liquid environments.The geometric criteria for hydrogen bond selections are the following: where θ is the OH • • • O angle and |d OO | is the distance between two oxygen atoms.A bi-exponential function has been fitted in HBACF to obtain the parameters corresponding to the short-time ( τ 1 ) and long-time ( τ 2 ) decay components of the autocorrelation function.
Finally, the lifetime of the hydrogen bond is calculated using Eq.(3): where τ HB is hydrogen bond lifetime.The last 20 ns trajectory was used for HBACF calculations. (1)

Results and discussion
It is well-known that confinement has a notable impact on water dynamics due to the change in HB structure of the molecules.In this work, we studied the effect of hydrophobic confinement on water dynamics inside different CNT sizes over a wide range of temperatures.We have computed the water radial distribution function (RDF) which is widely used, demonstrating the probable density of certain atoms as a function of distance.RDF is described by the probability of finding a particle at a distance r from a reference particle compared to the inhomogeneous distributions of the particle.The RDF for the oxygen atoms of the water molecules with respect to the CNT axis was calculated and shown in Fig. 2. The RDF has been computed using the following equation: where ρ is the number density.r i and r j represents the position of the i th and j th oxygen atoms of water molecule.
Figure 2 shows the RDF of confined water inside 1.0 nm, 2.0 nm, 3.0 nm, and 5.0 nm CNTs at different temperatures.Although the calculated results show slight variations among the different CNT sizes, the primary (6) www.nature.com/scientificreports/characteristic of the RDF curve remains consistent.The first peak corresponds to the first neighbor distance of the water molecule oxygen atom ∼ 2.84 Å.The amplitude of other peaks decreases rapidly with an increase in the distance between the water molecules and remains visible up to the third layer.1][22] .Furthermore, these results show that the ordering of water molecules is inversely proportional to temperature; the intensities of the RDF peaks decrease upon increasing the temperature.We find that the first peak position does not change with the increase in confinement size, however, the second peak shifts slightly with the increase in confinement size.The decrease in RDF peaks indicates a transition from the long-range structural order to a liquid-like phase.In addition, we also calculated the density map of oxygen atoms of water molecules in the XY plane.The density map was obtained by dividing the corresponding plane into the square bin of 0.1 Å length and then counting the number of oxygen atoms in each square bin.Higher oxygen densities are shown in red, while the low densities are shown in blue color.The inset of Fig. 2 shows the density map of water in each CNT at a temperature of 280 K.In all cases, density maps of water molecules are arranged in layered structures (rings in the density map), in agreement with the previous work 9,20,21 .The density map displays that upon increasing the CNT sizes the number of layers increases, consistent with RDF results.This layered structure observed in our simulated system does not depend on the temperature.It is clear now that confinement forces water molecules to form a layered structure but how the water ordering varies in these layered structures has not been discussed previously.
Recently, both with experimental results and MD simulation, we found that water molecules close to the center of CNT diffuse faster compared to the water molecules adjacent to the CNT wall, but the exact explanation of this observation remains elusive.To further address this, we have computed the local water order parameter.We relied on two widely used order parameters for liquids 37 i.e.Orientational tetrahedral and Translational tetrahedral order parameters.The orientational tetrahedral order parameter S q has been calculated for water inside different sizes of CNTs (Fig. 3A).As seen from Fig. 3A, S q in each CNT size is inversely proportional to temperature.This is expected since thermal energy provided to water molecules changes the structure from a quasi-icy structure toward a gas phase.At 300 K, the S q for the confined water in 1.0 nm CNT is 0.2744, 0.305 for 1.4 nm, 0.276 for 2 nm, 0.285 for 3 nm, and 0.359 for 5 nm.The results indicate that confined water molecules behave differently compared to the bulk water 37 .Results, excluding those in 2.0 nm indicate that S q is also inversely proportional to the CNT size.At first glance, this might be explained by the effect of the CNT wall on water molecules, at smaller CNTs compared to those in larger CNT sizes.It is important to point out that the analysis of Fig. 3A is based on considering the total number of water molecules inside CNTs.In Fig. 4, we observe the probability distribution of S q parameters versus temperature for all the CNT sizes considered in this study including the total number of water molecules.The distribution of this parameter for water molecules inside 1.0 nm and 2.0 nm shows two distribution peaks.The distribution in the case of 2.0 nm is shifted toward a low-value gas phase causing the overall weighted value of S q in 2.0 nm to give a smaller value compared to that of 1.0 nm as seen before (Fig. 3A).At larger CNT size the distribution gets broader and the overall value of S q decreases as also seen in Fig. 3A.
In Fig. 3B, the average translational tetrahedral order parameter ( S k ) of water molecules inside CNTs of various sizes is presented at different temperatures.Specifically, at 300 K, the S k for a 1.0 nm CNT is 0.905.As the CNT size increases, the S k values exhibit an upward trend: 0.988 for 2.0 nm, 0.992 for 3.0 nm, and 0.994 for 5.0 nm.We find that by increasing the confinement size, the S k value tends to approach the bulk value of 0.999, as reported in previous studies.This is because more water molecules fill the empty space, leading to a density approaching that of the bulk.Similar trends were observed in RDFs.Unlike S q , S k does not seem to offer a useful complementary measure of the change in the local structure of water molecules inside CNTs 37 .Since S k focuses on the variance of radial distances between an oxygen atom and the four nearest oxygen atoms, therefore, for 1.0 nm CNT, the S k value decreases with an increase in temperature.This is because at 260 K, the water molecules are ordered, radial distance fluctuations are less, and as the temperature is increased water molecules fluctuate more.Conversely, for the 2.0 nm CNT, the S k value decreases initially with temperature, but further temperature increments result in an increase.Notably, for 3.0 nm and 5.0 nm CNTs, the S k value shows minimal variation with temperature changes.
Previously, we reported 10,[20][21][22] , as well as other research groups 9 that water molecules inside CNT sizes rearrange themselves in coaxial cylindrical shapes, diffusing with different motilities along the axis of the CNTs.The self-diffusion coefficients versus temperature increase upon increasing CNT size, reaching maximum value when CNT size is about 2.0 nm -2.5 nm, then decrease and reach that of bulk phase at larger CNT sizes ( > 4.0 nm).The anomalous behavior of S q in the data of Fig. 3A at 2.0 nm might be related to water mobilization.To further investigate this, we have divided the volume inside each CNT channel into several regions, and then S q values www.nature.com/scientificreports/are calculated for each of those regions.For water molecules located within a cylindrical shell extending from the center of the CNT to 3.0 Å, results are presented in Fig. 5. Unlike the results presented in Fig. 3 (where the total number of water molecules was considered), the anomalous behavior observed inside 2.0 nm is no longer seen, since the CNT wall effect has been excluded.At each temperature, the S q is inversely proportional to the CNT size (including 2.0 nm data).The probability distribution S q of the water molecules from the center of the CNT to 3.0 Å (shown in Fig. 6) reflects this fact.Notice the difference between the results of 2.0 nm in this case compared to that of Fig. 4 where all the water molecules were considered for the calculation.This comparison confirms that inside 2.0 nm CNT size, water molecules close to the walls have a more structured quasi-icy phase compared to those at the center of the channels.This is further supported by observing the probability distribution of S q for water within larger CNT sizes.With ample space available, the structural arrangement of water molecules tends to lean towards a gas phase (lower values of S q ).To elucidate how the local ordering of water molecules affects the mobility of water molecules and hydrogen bonds, we have calculated both the self-diffusion coefficient ( D z ) and HBs at 300K.At this temperature, the D z value of the confined water in the 1 nm CNT is 0.882 × 10 −5 cm 2 /s .As the CNT size increases, the D z values increase 2.33 × 10 −5 cm 2 /s in 2nm CNT, 2.74 × 10 −5 cm 2 /s in 3 nm, and 2.42 × 10 −5 cm 2 /s in 5 nm, while the number of HB is 2.22 for water in 1 nm CNT, 2.61 in 2 nm, 2.71 i 3nm, and 2.80 in 5nm CNT.The result indicates that order parameters are directly proportional to the self-diffusion coefficient and HBs.
As the RDF results and density map show that water molecules inside the CNT channel of different diameters are arranged in a layered structure, it would be important to examine the dynamics of water in these layers.We have calculated the HBs for water molecules in each layers .Figure S1 shows the average number of hydrogen bonds per water molecule in various water layers inside differnet CNTs sizes.

Hydrogen bond dynamics
To evaluate the resilience of HBs among water molecules within the CNTs, we analyzed their longevity.This assessment included computing the HBACF, discussed in the method section, to investigate how both the size of confinement and temperature affect the stability of these crucial bonds.The total number of water molecules within the CNTs was considered.Figure 7 illustrates the variation of HBACF for different CNT sizes (1.0 nm, 2.0 nm, 3.0 nm, and 5.0 nm) at room temperatures.In Fig. 7's left panel, we observe an exponential decline in HBACF over simulation time across all CNT sizes.However, the rate of this decline noticeably varies with different confinement sizes.The inset within the figure distinctly highlights that water molecules confined within smaller CNTs (1.0 nm) sustain their HB networks for extended periods compared to those within larger CNTs.
The HB lifetime was determined by fitting the double exponential function as detailed in the method section.Figure 7B illustrates the HB lifetime across different CNT sizes at various temperatures, revealing a decreasing trend with rising temperature as well as CNT size.At 260K for smaller CNT (1nm), the HB lifetime is ∼ 3.3 indicating an ice-like structure which is consistent with our previous observation.In our previous study 22 , we calculated the HB lifetime using the TIP3P water model.At 300 K, the HB lifetime values for different CNT sizes were 1.0 ps (in 1.0 nm), 0.96 ps (in 2.0 nm), 0.94 ps (in 3.0 nm), and 0.91 ps (in 5.0 nm CNT).These values are Figure 5. Orientational tetrahedral order parameter ( S q ) of confined water molecules inside various CNT sizes, at different temperatures.To exclude the effect of the CNT walls, the analysis considers water molecules within a cylindrical shell ranging from 0 Å to 3 Å, with 0 Å referring to the center of CNT channel.Error bars represent the standard deviations of S q .notably smaller compared to the SPC/E water model.This shows that the HB network is less structured in the TIP3P model, emphasizing that the SPC/E water model is the optimal choice for the study of water dynamics.
As we know water molecules inside CNT are arranged in coaxial water cylinerical (CWC) sheets [20][21][22] .Therefore, we examine the HB dynamics in each CWC sheet using the HB lifetime.Table 1 presents the hydrogen bond information for water molecules corresponding to each CWC sheet inside different CNT sizes at temperatures 260 K, 280 K, 300 K, and 320 K.The table highlights the notable differences in water dynamics among the various CWC sheets.According to the data, at low temperatures, the outermost layer of water exhibits the longest HB lifetime providing a potential explanation for its slower diffusion, whereas as we increase the temperature, the HB lifetime remains the same throughout the different layers.These findings align with those of the previous studies.

Conclusion
In this study, we performed MD simulations to examine how both size and temperature influence the ordering of water and the dynamics of hydrogen bonds within confined CNTs.Our analysis spanned pore sizes from 1.0 nm to 5.0 nm across temperatures ranging from 260 K to 320 K.We focused on water molecule distribution and structural ordering of water in these confined environments.Additionally, we explored the impact of temperature Water molecules inside a cylindrical shell ranging from 0 Å to 3 Å are considered in this analysis.Note the difference between the results presented here and those in Fig. 4, where the total number of water molecules was considered.This distinction is particularly evident in the data for water confined within the 2.0 nm CNT. and confinement size on the hydrogen bonds and their lifetime.The radial distribution and density maps showed that water molecules forms cylindrical layers.To examine the ordering of watrer molecuels in these layers, we calculated the order parameters in various warter layers.The Orientational order parameter highlighted that in smaller nanotubes water adopts an ice-like structure near the tube walls, with the ordering pattern changing as the temperature rises.We find that water near the CNT center displayed a lower value of S q ( ∼ 0.34) compared to that of water mol- ecules near the wall ( ∼ 0.63).The higher S q value demonstrates that water molecules near the CNT wall form an ice-like structure.Due to this a significant reduction in the diffusion coefficient of water molecules close to the wall was observed both in experiments and MD simulations.Our findings for 2.0 nm CNT revealed the formation of an ice-like structure of water molecules adjacent to the CNT wall at low temperatures, a novel observation not previously reported.For 2.0 nm CNTs, the orientational order parameter ( S q ) exhibited anomalous behavior with increasing temperature, consistent with diffusion observations from prior studies.The hydrogen bond correlation function of water within CNTs exhibited a slower decay compared to bulk water, with the decay rate decreasing as CNT diameter increased.In larger CNTs, the hydrogen bond lifetime of the innermost layer was shorter than that of other layers and was temperature-dependent.We also demonstrate the water ordering CNT size above 3.0 nm is almost the same with that in the bulk.This study provides insights into the behavior of water molecules within a hydrophobic confinement.Overall, our study offers a detailed understanding of how confinement sizes and varying temperatures impact the structure and dynamic properties of water.These findings are important for applications in nanofluidics and other applications where water-CNT interactions are crucial.The discovery of ice-like structures and the anomalous temperature-dependent behavior of order parameters in specific CNT sizes provide new insights into the unique properties of confined water.

Figure 2 .
Figure 2. Radial distribution function of water Oxygen atoms inside various CNT sizes, at different temperatures (A) 1.0 nm, (B) 2.0 nm, (C) 3.0 nm and (D) 5.0 nm.The insets are the density map of water molecules in these CNTs, at 280 K, showing the formation of the layered structure.The high-water-density regions are shown in red color, and the low-density regions are shown in dark blue color.

Figure 3 .
Figure 3. (A) Orientational tetrahedral order parameter ( S q ) of confined water molecules inside various CNT sizes, at different temperatures, (B) Translational tetrahedral order parameter ( S k ) of confined water molecules inside various CNT sizes, at different temperatures.S q and S k are obtained by considering the total number of water molecules within the CNT channels.Error bars indicate the standard deviations of the both parameters.

Figure 4 .
Figure 4. Probability distribution of Orientational tetrahedral order parameter of confined water molecules inside various CNT sizes, at different temperatures (A) 1.0 nm, (B) 2.0 nm, (C) 3.0 nm and (D) 5.0 nm.The analysis takes into account the total number of water molecules confined within the CNT channel.

Figure 6 .
Figure 6.Probability distribution of Orientational tetrahedral order parameter of confined water molecules inside various CNT sizes, at different temperatures (A) 1.0 nm, (B) 2.0 nm, (C) 3.0 nm and (D) 5.0 nm.Water molecules inside a cylindrical shell ranging from 0 Å to 3 Å are considered in this analysis.Note the difference between the results presented here and those in Fig.4, where the total number of water molecules was considered.This distinction is particularly evident in the data for water confined within the 2.0 nm CNT.

Figure 7 .
Figure 7. (A) Hydrogen bond auto-correlation function of confined water molecules inside various CNT sizes, 300 K.The inset is the magnification of the data between 0 to 5 ps.(B) Temperature dependence of the hydrogen bond lifetime of the confined water molecules, at different temperatures.

Table 1 .
The hydrogen bond lifetime of water molecules within different CNT sizes, at different temperatures.The Water shells in the table correspond to the coaxial water cylindrical sheets formed in CNTs due to water confinement.